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I. INTRODUCTION 

There are essentially two reasons that can cause a cur- 
rent of particles like, e.g., electrons through solids: Either 
there is an electric field dragging the electrons in some 
direction or the electron density is spatially non-uniform, 
i.e., there is a density gradient causing the electrons to 
diffuse. For systems featuring normal transport the cur- 
rents in those two cases are supposed to be determined 
by 

j = L F F , (la) 



j = -L D Vp, (lb) 

where j is the current, F the external force, p the spa- 
tial density and the L's are the pertinent response coeffi- 
cients. As will be briefly outlined below the derivation of 
(|lafl from the corresponding scenario is rather straight- 
forward and leads to the Kubo formula^ (KF). A di- 
rect derivation of i|lb|l from its underlying scenario seems 
to be significantly more subtle. This is a crucial point 
for the analysis of heat conduction since, although being 
rather similar to the above case otherwise, the scenario 
corresponding to (|laf> does not exist for heat transport. 
There simply is no external field which could exert a force 
on heat, heat is always driven by an energy density (tem- 
perature) gradient as described by ljlb|) . Nevertheless, as 
will also be sketched below, there have been attempts to 
apply the derivation of l|la|) to thermal conduction which 
eventually implies the application of the KF also in this 
cas o^Vi 4 . 

Let us briefly recall the derivation of the KF. An exter- 
nal field gives rise to an additional energy H' = J pUdV 
with — VC/ = F . This term is routinely treated as a per- 
turbing addend to the Hamiltonian of the system yielding 
an expression for the induced current. In this expression 



H' itself no longer appears only its derivative with re- 
spect to time which in the case of a (spatially constant) 
external field reads H 1 = — jF. Thus one obtains a rela- 
tion of the form lllal^ . Lp is given by the Kubo formula 
and reads 

L{u) = - die- 1 "* drTr{p j(0)j(t + ir)}, (2) 
v Jo Jo 

where po is the Gibbsian equilibrium state, V the volume 
and L(u>) is the response coefficient which describes the 
conductivity at frequency ui. 

In the case of heat conduction no additional energy 
arises from the internal force (temperature gradient) and 
thus the stimulus of the current cannot be incorporated 
into the Hamiltonian. Hence the above line of reason does 
not apply. To nevertheless justify a Kubo-type formula 
for thermal conductivity basically two types of arguments 
are brought forth: 

i) The hypothesis of local equilibrium 
The state that one gets by boldly writing down a Gibbs 
state with a spatially non-uniform temperature that 
varies little (AT(x)) around some mean [3q — 1/Tq reads 

p loq = ^ 1 exp( -pojdx 1 - ^ T(3;) /i(x)) (3) 

(Z being the partition function, h(x) the local energy 
density and T respectively (3$ some mean temperature) 
and is called a local equilibrium state (cf. Ref •£*&). Its 
physical significance is somewhat vague since it is not 
a real equilibrium state. However, adding a "pseudo- 
perturbation" of the form H' ps = - J dxAT(x)h(x)/T 
to the system's Hamiltonian and calculating the Gibbs 
state at uniform (3o of the perturbed system formally 
yields the above local equilibrium state (0. Thus this 
term is said to somehow model the effect of the inter- 
nal temperature gradient, although the local equilibrium 
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state features no current. Nevertheless, proceeding like 
described above and taking H' ps for an external pertur- 
bation yields a transport coefficient as given by J5J) (just 
multiplied by (3q) with a thermal current j defined by 
h + Vj = 0. 

ii) The entropy production argument 
If, e.g., an electrical current in some conducting solid 
runs along an electric field, potential energy arising from 
this field, H' , is converted to heat, Q. Thus one has 
—H' = Q. Hence the entropy production is S > Q/T. 
On the other hand entropy production is assumed to be 
of the form S = J ■ VT/T 2 (see Refsi7i 8 i 9 i ln ). Combin- 
ing these equations and boldly replacing the ">" by an 
"=" yields H' = -J-VT/T. This can formally be incor- 
porated in the derivation of the KF outlined above and 
yields the same result as the hypothesis of local equilib- 
rium. The crucial shortcoming of this argument is that 
if heat flows along a temperature gradient, obviously no 
heat is produced, only entropy. Thus, in this case one 
would definitely have to consider a ">" and then this 
argument yields no concise result. 

Thus it has often been pointed out that those con- 
cepts do not provide a rigorous justification of the KF 
for thermal conductior^*ii*i2iiiiiiSii& which remains an 
open question. Furthermore, the KF has been counter- 
checked only for few concrete systems, see RefsiiZ*i&. De- 
spite those conceptual problems the application of the KF 
to heat conduction is today a standard techniqusiSiSSiSi. 
Thus the main intend of the paper at hand is to give a 
derivation of a Kubo-type formula for thermal conduc- 
tion which is more convincing, thereby also pointing out 
the limits of its validity. 



II. KUBO FORMULA FOR FINITE QUANTUM 
SYSTEMS 

Before we proceed with our main line of thought in 
Sect. IIIII we shortly comment on another difficulty that 
arises if one wants to compute the transport coefficients 
of finite quantum systems from the KF. This point de- 
serves consideration, since, in general, an infinite sys- 
tem cannot be treated numerically, thus in many cases a 
transport coefficient of an (periodic) infinite system can 
only be computed from analyzing a finite piece of the 
infinite system (cf. Refi22). 

Of course, a key question is whether or not the trans- 
port is normal at all. For infinite systems the coefficient 
is believed to take the form L(lu) = Dd r S(0) + k{lo). Dd r 
is called the Drude weight and whenever it is nonzero for 
infinite systems the transport is assumed to be ballistic. 
If it is zero, the normal (diffusive) conductivity is sup- 
posed to be given by k(0), k(u>) being a smooth function 
without any singularitiesSS^. However, this distinction 
is problematic if the KF is evaluated on the basis of a 
finite quantum system. In this case L(u>) consists of a 
sum of delta peaks at different frequencies without any 
non-singular contribution22*2£. So, technically speaking, 



one cannot find normal transport in a finite quantum 
system, it is either ballistic or none. There are differ- 
ent ideas on how the transition to infinite systems could 
produce normal transport. Sometimes the singularity at 
lu = is assumed to broaden ( "imaginary broadening" ) 
but maintain its weight such that L(0) = Dd r /T is non- 
singular, where r is some inverse width of the Drude-peak 
which is hard to evaluate from a closed finite system^. In 
other approaches L(u>) is averaged over small frequency 
intervals Slu and the resulting smooth function is extrap- 
olated down to lu = in order to determine the normal 
conductivity22i2ii. The result depends, of course, on how 
5lu is chosen. In general it appears to be difficult to de- 
termine transport type and conductivity merely from the 
Hamiltonian of a finite piece of an infinite system without 
any further assumptions. 

Thus, except for demonstrating the limited validity of 
the KF, the paper at hand also suggests a consistent 
method to infer the heat conduction behavior of an in- 
finite system from analyzing an adequate finite piece of 
it. 



III. MODULAR QUANTUM SYSTEMS AND 
DIFFUSIVE DYNAMICS 

We consider systems consisting of identical many-level 
subunits which are weakly coupled by identical next- 
neighbor interactions. We call those systems "modular" . 
For simplicity, we analyze chains and rings of that kind. 
Their Hamiltonians may be denoted as 

N N-l 

ff=^/i(/i) + ^t>(/i,M + l), (4) 
fi— l /^— l 

where h(fj.) is the local Hamiltonian of some subunit /i, N 
the total number of subunits and V(ft, [/, + 1) represents a 
next-neighbor interaction between the subsystems [x and 
/i + f (Schematic examples of such modular systems are 
depicted in Fig.^and Fig. EI) • The one-dimensional char- 
acter is not crucial here. Everything derived below can 
be generalized straightforwardly to modular systems fea- 
turing arbitrary multi-dimensional "net-structures" . So 
what sort of physical systems are modular systems? First 
of all interacting nano-structures like arrays of quantum 
dots, etc., might fit this scheme. But modular structures 
may also be achieved by operationally coarse graining 
periodic (or also slightly disordered) systems like spin 
chains, crystals, etc. in modules such that each module 
containins many elementary cells. The interactions be- 
tween the modules then represent the couplings. Since 
interactions are typically rather short ranged, increasing 
the "grain size" will eventually result in a description 
in which only adjacent modules are coupled. Further- 
more, these next neighbor couplings will become weaker, 
such that they finally might be considered as weak. This 
"weak coupling" , as well as other criteria which deter- 
mine whether or not there is regular transport, will be 
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given more precisely below. They may all depend on the 
choice of the grain size. This reflects the fact that regu- 
lar thermodynamic behavior may only be expected for a 
spatially coarse enough description. Whether or not the 
criteria for regular transport are fulfilled at some grain 
size has to be investigated for any system individually. 
However, the description of a system as a modular sys- 
tem should be possible for very many systems. Thus, 
in order to keep everything as general as possible, we 
do not specify our systems in much more detail than al- 
ready given by |@J (not even the numerical examples in 
Sect. lVlill . But the concrete application of the results at 
hand to, e.g., spin chains is under way. 

What does diffusive transport mean in the context of 
modular systems? Diffusive heat transport is defined by 
Fourier's law which is essentially given by (|lb|l . For mod- 
ular systems we replace the spatial gradient by the en- 
ergy difference between two adjacent subunits, i.e., the 
current from subunit fi to subunit /i+ 1 should obey 
= D(h([i) — h([i+l)). A discrete form of the continu- 
ity eq. for modular systems reads /i(/i) = j'Qi — 1) — 
Combining those eqs. yields, e.g., for a chain 

^±=D(h(2)-h(l)), 

^hL = D(h((i-l)+h(n + l)-2h( t i)), (5) 
^L = D(h(N-l)-h(N)). 

(This may be viewed as a discrete form of the time- 
dependent version of Fourier's law: h(x) = (k I 'c) Ah(x) .) 
If the motion of energy through the closed system 
(which is entirely controlled by its Hamiltonian and the 
Schrodinger eq.) can be described by the above set of 
eqs., the thermal transport is diffusive. The conductivity 
is obviously related to D. 

Our paper is roughly organized as follows: We consider 
the dynamics of a modular system without any external 
forces, starting from a local equilibrium state as given 
by @. We find that, under various conditions on the 
model parameters, the motion of the energy can be de- 
scribed by JSJ for a short first time-step (cf. Sect. IV}l . 
where D is essentially determined by the KF. After that 
time-step the system is unfortunately no longer in a local 
equilibrium state. It can, however, be shown that almost 
all states sharing some crucial properties with the local 
equilibrium state will also give rise to energy dynamics in 
accord with JSJ for another time-step. That means, fully 
evolved local equilibrium is dispensable. Thus one can, 
iterating the result for short time-steps, conclude that © 
provides a correct description for all times (cf. Sect. IVl|) 
and thus justify the application of the KF. Eventually 
we evaluate the KF for some concrete finite models and 
compare the results with the energy dynamics obtained 
from a direct numerical integration of the corresponding 
time-dependent Schrodinger eqs. (cf. Sect. IVlIj l. 



IV. DEFINITION OF LOCAL ENERGY 
CURRENTS IN MODULAR QUANTUM 
SYSTEMS 

Consider a Hamiltonian as given by Q. In the paper 
at hand we define the local energy operator at site fi sim- 
ply as h(p) rather than h(/j,)+(V(fi, [i—l)+V(fi, fi+l))/2. 
This way the sum of all local energies does not represent 
the total energy and is not a strictly conserved quantity. 
But the local energy operators are defined on strictly sep- 
arated subspaces of the product space on which the full 
system is defined. Excluding the interactions from the 
local energies makes a consistent partition of the system 
into mutually disjoint (smallest) subunits possible, such 
that on each subunit a local (equilibrium) state may be 
defined independently. This would be impossible if one 
included the interactions in the local energies. However, 
if one wants to define an energy current based on the 
evolutions of those local energies, the sum of the local 
energies has to be at least approximately conserved, i.e., 
the part of the full energy associated with the local en- 
ergies, Tv{ph(fi)}, has to be much bigger than the part 
associated with the interactions, Tr{/5 V(fi, fj, + 1)}. This 
eventually means that the interactions have to be weak. 

For conserved quantities the discrete continuity equa- 
tion h(n) = j(fx — 1) — suggests the definition of a 
local current operatofli5iMii22 via (h = 1) 

±h(jM)=i[H,hM] (6) 
= i[V( t i - 1, fj)Mri] + i\V{l*, M + l),h(j*)] . 

S * ' S v ' 

However, rewriting the above equation for the temporal 
change of h(^i + 1) produces another expression for the 
local current j(fx) — i[V(fi,(x + l),h(fx + 1)]. This can 
only be consistent if 

\y{n,H + l),hQi) + h(n + l)]K>Q. (7) 

But this expression can be interpreted as the tempo- 
ral change of the sum of two adjacent local energies, if 
only the interaction between the respective subunits was 
present. According to the weak coupling precondition, 
this sum of local energies is approximately conserved. Or 
at least its fluctuations are small compared to its mean 
value as long as the energy contained in the subunits is 
not too small. This means in particular that Q may 
safely be assumed as long as the temperature is not too 
small (we will come back to that condition later). Thus 
a sort of symmetrized local current is defined by 

j(/i)=i[t>(/i,/i+l),S], §=±(h(p + l)-hW). (8) 

The local current operator is strictly defined within the 
product space spanned by the corresponding adjacent 
subunits. Its expectation value is determined only by the 
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reduced state of those adjacent subunits. Thus, within 
this framework, the relation between current and tem- 
perature difference may be determined on the basis of a 
reduced system consisting of only two interacting sub- 
units. In this reduced system the current is simply the 
temporal change of S. And S may be interpreted as the 
operator for the skewness of the energy distribution be- 
tween the two subunits. 



V. LOCAL ENERGY CURRENTS AS A 
RESPONSE TO LOCAL TEMPERATURE 
GRADIENTS 

We now investigate the relation between temperature 
difference AT and the short time behavior of the energy 
current j. Inspite of having commented in Sect. [I] on 
the local equilibrium state in a rather critical way we 
now analyze the short time dynamics, i.e., the forma- 
tion of a current of an initial local equilibrium state as 
given by J2J. Since we are not deducing any pseudo po- 
tentials here, this does not imply the application of the 
hypothesis of local equilibrium in the sense described in 
Sect. [I] However, if one starts with a local equilibrium 
state and considers only a short time step, the question 
arises whether the system will be in (another) local equi- 
librium state after this time step? The answer is no, but 
this issue is addressed in detail in Sect. IVII For the mo- 
ment we thus simply consider the local equilibrium state 
p(T, AT) which is defined for the two coupled subunits 
which form our reduced system (for simplicity, those are 
only labeled "1" and "2" in the following): 

-(fe(l)+fe(2)) 
T 



* AT - ex P 
p(T, AT) w p (l + —S) , po := 

(9) 

[Z is the partition function of one subsystem and we 
chose units of temperature and energy as fee = 1, 
h = 1.) Here po is obviously a "global" equilibrium 
state with both subunits at the same temperature T. 
Unfortunately, the current for this state vanishes, i.e., 
Tr{p(T, AT) j} = (if the partial traces of V with re- 
spect to the subunits vanish, which can be demanded 
without loss of generality). 

Thus one has to proceed in a slightly different way: 
We analyze an approximate short time evolution of the 
expectation value 



S(t) = Tt{p(T,AT)S(t)} 



(10) 



and consider it's derivative with respect to time for some 
small but finite r. This derivative is, according to the 
Heisenberg equation of motion and the definition in JSJ , 
the current at time r. The time evolution of the operator 
S(t) is computed by means of a truncated Dyson series, 
i.e., 



S(t) = D\t)SD(t) with 
D(t) ^ 1 - iU^r) - U 2 (t) 



and the time evolution operators (see KeiM^) 

#i(r) = f T dr'V(r'), (12) 



U 2 (r) = \ dr' [ T dr" V(r') V(r") . (13) 
Jo Jo 

Here the time dependence of the V(r)'s is defined on the 
basis of an interaction picture, i.e., only generated by 
the local Hamiltonians h(l) + h(2) rather than the full 
Hamiltonian. This yields in second order approximation 



S(t) « S+ ipTi, S] + UxSUt - (SU 2 + U\S) . (14) 

(For simplicity of notation, we suppress the time argu- 
ment of the U's here and in the following.) Computing 
S(t) defined in l|10(l by plugging in the operator (|I4|I and 
the state © yields 



S(r) 



^(tt{poS 2 } + ^Tt{p [Ui,S} 2 }) . (I') 



where we have exploited U 2 + U 2 = U 2 (which follows 
from the definitions HI2|I and C3J) and [pojS 1 ] = as 
well as [E/i,/3o] ~ 0. As already mentioned the latter is 
valid for high enough temperatures (remember discussion 
below Q). Realizing, by using the definition of U\ and 
(H, that 



[U u S] = -i / j(r')dr' 
Jo 



(16) 



where again j(r') is defined according to the interaction 
picture, one can write the derivative with respect to time 
of (US!) as 

±S(t) = -^J\r{p j(r')j(r)}dr' . (17) 
Eventually, substituting t' = r — r' and t = t 

±s(t)=m = J \T{p j(o))(t')} dt' , as) 

= :C(t') 

the current is essentially given by an integration over the 
current auto-correlation function, C(t'), just like in the 
KF. Let the timescale on which this correlation typically 
decays be r c . Then, if the approximation for S(t) (sec- 
ond order Dyson scries) holds for times larger than r c 
(which will be analyzed below), the current will indeed 
assume a steady value after r c which is proportional to 
AT. Thus, first of all, in the case of free heat transport 
without any external baths Fourier's law is confirmed for 
the short-time evolution of a local equilibrium state. The 
conductivity is now defined by 



(11) 



i r° 

K = f2 J C(t)dt with t c < t < T d 



(19) 
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where Td is the timescale on which the second order ap- 
proximation for the Dyson series brakes down (r<j will be 
evaluated more concretely below). Within this interval 
the integral should not depend much on its upper limit 
tq. Nevertheless, for finite quantum systems (for which 
C(oj) is just a set of peaks) this value may differ consider- 
ably from the one which is produced by letting the upper 
limit go to infinity. To see this and the relation to the 
KF it is instructive to consider C(oj) and C(u>) defined 
by 



C{w) 



2tt 



■ C(t)dt 



CV) = 1- / C(w)6w 



(20) 
(21) 



Obviously, C(ui) is the Fourier transform of C(t) and 
C(uj) a slightly "smeared out" version of it. Let C{t) be 
the back transform of C{lo). Now, as long as t <C 2tt/Su>, 
one has C(t) « C(t). Thus 



To pTQ 

C(t)diw / C(t)dt 
Jo 



(22) 



as long as To <C 2n fbio. Since C{t) (other than C(t)) does 
not feature any recurrence time if C (w) is a smooth func- 
tion of frequency, we may now drive the upper boundary 
of the second integral to infinity without much changing 
its value. This yields 



C(t)dt : 



C(t)dt = ttC(lu = 0) 



(23) 



(In the following we denote C(u> = 0) simply as C(0).) 
This, however, can only hold if r c <C 2tt/Suj. A rough 
estimation for r c is given by r c rj 2ir/ Aw where Aw is the 
width of the spectrum of C(u>). Consequently, Aw 3> <5w 
has to hold. Of course, Sui can always be chosen to fulfill 
this, but if it is chosen too small, C(ui) is not necessary a 
smooth function of frequency. Thus 123|) eventually holds 
if C'(w) has a reasonably well defined "peak density" on 
a frequency scale small compared to its width Aw. In 
this case C(0) does not depend much on 5u> and there is 
no need to define 5u> with extreme precision. 

The KF (also evaluated for a reduced two-subunit sys- 
tem with V, as routinely done, replaced by the number of 
contacts between identical subunits, i.e., V = 1) and its 
limit for w <C T, L'(u), read in terms of the correlation 
function 

L(u) = n 1 ~ °~J /T C(u) , L\w) = ^C{u>). (24) 

This is to be compared with the conductivity as obtained 
from (UU, J2U), fgft which reads 



(25) 



Obviously, n formally equals L(0), the peak density of 
L(lu) at frequency zero. As mentioned in Sect. [H] this 



quantity has been suggested to compute the conductiv- 
ity of infinite systems from finite models. In this sense 
and for t < Td the KF is valid for the computation of 
thermal conductivity, although there is no external po- 
tential in this scenario. Note that the correlation func- 
tion in the case of the KF is given by the full Heisenberg 
dynamics of the system based on the full Hamiltonian, 
whereas n is determined by the correlation function based 
on the interaction picture dynamics, i.e., without taking 
the interaction into account. However, since the whole 
theory is formulated for the weak coupling limit, this is 
not going to make a big difference as will be numerically 
demonstrated below (cf. Sect. IVIlJl . 

So far we have always assumed that r c < Td- But is 
that justified? The first order term of the Dyson series is 
iU\. Thus an estimate for the magnitude of the first order 
term which puts weight to the energy subspaces being 
proportional to their occupation probability is given by 



F(t) := Tr{p U?} 
Exploiting 1)12(1 one finds 



(26) 



T F(t)tt / Tr{p o V(0)V(t')}dt' . (27) 
dt J . ' 

— Cy(t') 

Since the decay of the interaction auto-correlation CV(i') 
will roughly proceed on the same timescale as the current 
auto-correlation, the temporal change of F(t) will assume 
a steady value after r c . Thus, after r c , one has F(t) ~ 
Cv(0)t, where Cy(0) is defined analogous to (7(0) but 
based on the interaction auto-correlation. The time for 
a valid approximation according to the Dyson series may 
now be defined as 



1 



Td 



C v (0) ' 



(28) 



a time for which F(rd) — 1 and consequently the second 
order approximation definitely brakes down. 

If, based on this definition, r c < Td is not fulfilled 
the current cannot be shown to assume a steady value 
proportional to AT and thus Fourier's law will, in gen- 
eral, not be fulfilled. That means the model does not 
show normal transport. According to (|27|l this happens 
if the coupling becomes to strong. This fact might be 
missed by simply evaluating the KF: Since C(u>) com- 
pletely scales with the square of the overall interaction 
strength (at least for C(ui) evaluated in the interaction 
picture) , the distinction between normal and non- normal 
transport cannot be made by simply looking at the fea- 
tures of C(u>) (this will be demonstrated in more detail 
in Sect. IVII|) . The transport behavior can also become 
non-normal if the coupling becomes too weak. Namely, 
for Td 3> 2it/5lo the current is well predicted by the trun- 
cated Dyson dynamics for times for which C(t) w C(t) 
does not hold any longer. This also means that the cur- 
rent cannot be predicted in the described way and thus 
Fourier's law will typically not be fulfilled (this will also 
be demonstrated in more detail in Sect. rvTT|) . 
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VI. HILBERT SPACE AVERAGE METHOD 
AND ITERATION SCHEME 

So far we have shown that, under given conditions on 
the model, for a short time period an energy current will 
flow between the subsystems which is proportional to the 
local temperature difference AT. But this has been en- 
tirely derived under the assumption that the initial state 
was a local equilibrium state of the type p(T, AT). So 
what happens after r^? It is tempting to iterate the pro- 
cedure, assuming that the state after n would be again a 
local equilibrium state only with some reduced AT. Un- 
fortunately, this is rigorously not the case. At t = the 
initial state p(T, AT) factorizes. This will for t > t<j 
no longer be the case. The state at t > features 
an instantaneous current, p(T, AT) does not. The full 
system's von-Neumann entropy always equals its initial 
entropy, but a local equilibrium state with reduced AT 
would feature a higher entropy. 

Thus, in the following we abandon the local equilib- 
rium state. Instead we show, that essentially almost all 
states that feature a certain S induce, for a short time- 
step, a current that is proportional to S. This, of course, 
results in a state featuring an accordingly reduced S. As- 
suming that this state also belongs to the above class, the 
above short time-step dynamics for S may be iterated. 
This yields continuous dynamics for S and hence for the 
current. The assumption is only reliable if the above class 
contains basically all states in quest. 

The formalism to implement this scheme is called 
Hilbert space average method (HAM). The key idea of 
this method is to replace the expectation values for ob- 
servables A of actual pure states (i{)\A\ip) (in our case 
the expectation value of the energy skewness (ij;\S\ip)) 
by their Hilbert space averages 

(29) 

This expression stands for the average of (ip\A\ip) over 
all that feature (ip\B\i/}) — b but are uniformly dis- 
tributed otherwise. Uniformly distributed means invari- 
ant with respect to all unitary transformations that leave 
(ip\B\i/j) = b unchanged. The replacement of actual ex- 
pectation values by their Hilbert space averages is only a 
justified guess if almost all individual \ip) yield expecta- 
tion values close to the Hilbert space average of the ob- 
servable. It can be shown that this is the case if the spec- 
tral width of A is not too large and A is high-dimensional. 
Full explanation of HAM is beyond the scope of this text 
and can be found in KeiM^. Here HAM is only to be 
applied. 

For the moment let the system be in a pure state. 
Then the current is given by the temporal change of 
(if)\S(r)\i/>}. Assume that the following set of expecta- 
tion values is known: (^?|P(?7)|^) = P(f]), where P(r)) 
is a projector, projecting out the subspace which corre- 
sponds to an energy interval of width Ar/ around E = r), 



i.e., 

P (V)= E I>' S X S '*I' ( 3 °) 

E=7]-Ar)/2 s 

where E are eigenvalues of E := h(l) + h(2) (sum of 
local energies) and s eigenvalues of S (energy skewness). 
Since, for high enough temperatures, the sum of the local 
energies is an approximately conserved quantity, the P{rf) 
are, for large enough A77, also approximately conserved. 
Assume furthermore that the initial energy skewnesses 
within all subspaces 77, i.e., {tp\P(ri)S\ip) — S(r)) are also 
known. Without taking any further information on \tp) 
into account the best guess on the evolution of (iP\S(t)\i/j) 
is given by the Hilbert space average 

li' l P\^( T )\' l P)i{(^\p(ri)\^)=P(r l ),{i,\P(ri)SW=S(r 1 )} = Tr{Sa} 

(31) 

with 

a = [IV')(V ; l]{(^|p( ?? )|^) = p(^) ! ^|p(^)5|^) = 5(^)} • (32) 

Now what is the above Hilbert space average a? Any 
unitary transformation G that leaves P(r]) and P(r))S 
invariant has to leave a invariant, i.e., 

e i6 de~ i6 = q with 

[G,P(r 1 )S} = [G,P(r 1 )} = 0. (33) 

This, however, can only be fulfilled if [G, a] = for any 
G. Furthermore, since the Hilbert space average a is to 
be computed under some restrictions (see (I32H '|. one has 
the following conditions 

Tv{aP( v )} = P( v ) , Tr{aP(r,)S} = S(r,) . (34) 

According to the invariance properties of a and the prop- 
erties of the operators P(rf) and P(rf)S (I33f) . one may thus 
write a as 

^\ Wy Tr{T(r,)} Tr{P(t))S*} J 

Defining p (r]) := P(?7)P(r;)/Tr{P(77)}, (33 may be 
rewritten as 

By choosing P{rj) and S{rj) to be equal to the expectation 
values of p(T, AT) for the respective operators, i.e., 

P( V ) := Tr{P(ri)p(T, AT)} = Tr{^P(v)} (37) 

AT 

S(r,) := Tv{P(r])Sp(T, AT)} = —Tv{p (v)S 2 } 

one gets J2 V Po(jl) ~ Po and thus a p(T, AT) as can be 
seen from comparison with Hence the Hilbert space 
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average over all pure states featuring the same expecta- 
tion values for P(rj) and P(r})S as p(T, AT) is p(T, AT) 
itself. Therefore 



(ij\S(TM) ^Tr{p(T,AT)S(r)} , 



(38) 



i.e., any pure state (regardless of whether it is entangled 
with respect to the subunits or whether its local entropy 
is maximum) featuring the same expectation values for 
P(rj) and P(r])S as the local equilibrium state p(T, AT) 
is most likely to yield the same local current as the lo- 
cal equilibrium state. Since incoherent mixtures of pure 
states featuring the same expectation values for P(rf) and 
P(i])S as the local equilibrium state are even closer to 
the latter then pure states, they are even more likely to 
induce the same current as the local equilibrium state. 
This holds for the set of P(rj), S(rj) as given by l(37j). But 
how will this set look like after tq, i.e., how do the P{r}) 
and S(ri) evolve during the time tq? Since the P(j]) are 
approximate constants of motion, they remain invariant. 
For the S(rj) we find from applying the scheme developed 
in Sect. E| to Tr {&P(t])S(t)} [cf. JTBJ, {El, 123, (HZ)] 



S(r),t + T ) aSfat) 



7rC(»7,0) 
Tr{Mv)S 2 } 

= :2D(rt) 



S(r,,t)T , (39) 



where C(r], 0) is analogous to (7(0) as defined by and 
(|18|1 but po replaced by po(rf). If t ne criteria from Sect. El 
are fulfilled and tq is comparatively short, this may now 
be iterated yielding 



df 



S( v ) = -2D(r,)S(r)) . 



(40) 



Thus, in general, any energy subspace has its own en- 
ergy diffusion coefficient. Therefore it might, be impos- 
sible to describe the energy diffusion behavior between 
the subunits entirely by one overall conductivity as im- 
plied by the KF, even if all the criteria from l|18|l are ful- 
filled. However, if the D(r]) corresponding to the energy 
subspaces that are occupied with significant weight and 
contribute significantly to the transport feature similar 
values, i.e, D(rj) w D, one gets only one single diffusion 
coefficient that then reads 



D 



nC(0) 
Tr{/5 5 2 } 



(41) 



(Note that Y^^C^q.Q) = In this case one might 

sum (|40ll over 77 which yields 



— S= -2DS 
dt 



(42) 



It is straightforward to show that the heat capacity c for 
one subsystem is given by c = 2Tr{poS 2 }/T 2 . Exploiting 
this and l|41|l one may rewrite Q42|l as 

^=~* S ° r 3{l) = - c {h{l)-h{2)) (43) 



1 
Se 

T 



AE 



® 



p = l 



p = N 



FIG. 1: Simple model to analyze diffusive energy transport: 
TV coupled subunits featuring a single ground state and an 
excitation band of n equally distributed levels. 



which, generalized to N subsystems, identifying D = k/c 
and combined with the discrete continuity eq. , yields JSJl . 
Hence this closes the loop for a microscopic derivation of 
(|5|). This essentially means that if the criteria for To, r c , Td 
from Sect. are fulfilled and the relevant D(rf) are sim- 
ilar, energy will diffuse from subunit to subunit as pre- 
dicted by Fourier's law. 



VII. APPLICATION TO MODELS 

In this Section the previously derived theoretical re- 
sults are compared with numerical data in the follow- 
ing way: For concretely defined models the energy dif- 
fusion coefficient D — k/c which appears in Q is com- 
puted from the KF. Then is solved for some non- 
uniform initial energy distribution. Furthermore, the 
time-dependent Schrodinger eq. is solved numerically for 
an initial state corresponding to the above initial energy 
distribution. From those data the exact evolution of the 
local energies h(p) = (ip(t)\h(p)\ip(t)) is computed and 
compared with the above solution of Only if there 
is good agreement, the system exhibits normal transport 
that may be characterized by the conductivity obtained 
from the KF. 

We introduce two classes of models which are primarily 
designed to represent modular systems in general rather 
than real physical systems. (For the impact on real phys- 
ical systems see Sect. Hill ) The first class of models fea- 
tures subunits with non-degenerate ground states, large 
energy gaps (AE) and one, comparatively narrow energy 
band (Se) each, as depicted in Fig. ^ Within one band 
there are n states featuring equidistant level spacing. The 
next neighbor interactions are defined as 

V = A «(*, ])P + (h l*)P- (J-, M + 1) + h.c (44) 

(h.c. stands for the hermitian conjugate of the previous 
sum.) Here P + (i,p) corresponds to a transition of the 
/i'th subunit from its ground state to the i'th eigenstate of 
the band. P~ corresponds to the respective downwards 
transition. The v(i,j) are randomly distributed complex 
numbers normalized to Y]^ ^ \v(i,j)\ 2 /n 2 = 1. Thus A 
sets the overall interaction strength. Due to this model 
design only one energy subspace contributes to transport 
at all, which is the "one-excitation subspace" defined by 
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FIG. 2: Heat conductivity over frequency for a ringlike system 
as depicted in Fig. (N = 6, other parameters see text). 
Crosses refer to the coupled system, as intended by the Kubo 
formula, and x's to the decoupled system, as suggested in 
Sect. El 



an overall energy E with AE < E < AE + Se. The 
dimension M = nN of this relevant subspace grows lin- 
early rather than exponentially with the number of sub- 
units. This allows to numerically analyze models with 
high enough n to fulfill all criteria for diffusive trans- 
port, but also up to fifteen subunits. Thus those models 
are meant to check whether or not the locally computed 
conductivity holds for arbitrarily many subunits, as pre- 
dicted by the theory. (For a "stand alone" treatment of 
such systems, cf. Refill.) 

We consider a ring of N — 6 subsystems with n = 500, 
A = 5 ■ 10~ 5 , AE = 10 and Se = 0.05. We find that with 
a frequency averaging interval of Slu m 10~ 3 = <5e/50 all 
the conditions mentioned in Sect. are fulfilled. Since 
only one energy subspace contributes to transport, none 
of the difficulties concerning different L(rf) discussed in 
Sect. I VII arises. Thus the KF may be evaluated at 
any temperature to compute D = k/c. We evaluated 
k(w) = C{u)/T 2 (cf. Sect.© for T = 1.4 on the basis 
of the coupled system, as intended by the Kubo formula, 
and on the basis of the uncoupled system, as intended by 
our argument in Sect. ]V\ Both results are displayed in 
Fig.0 Obviously, there is a good agreement between the 
two graphs, it appears to be irrelevant whether the cor- 
relation function is evaluated on the basis of the coupled 
or the decoupled system. (It should be mentioned that 
in both cases C(uj) features finite contributions exactly 
at uj = 0, thus there is a finite Drude peak.) From Fig. [2 
we find k = 1.6- 10~ 3 (oj = 0) and calculating the specific 
heat for one subunit yields c(T = 1.4) = 10.5. This yields 
D = 3.142-10~ 4 . The corresponding solution of (J5J for all 
energy initially concentrated in one subsystem is shown 
in Fig. [2 Furthermore, we solved the Schrodinger eq. for 
a corresponding pure initial product state, featuring one 
subunit in a randomly generated state restricted to the 
excitation band, and all other subunits in their ground 
states. The result is also shown in Fig. [3] Obviously, 
there is fairly good agreement. We checked rings up to 
fifteen subunits and always found good agreement. So far 
the KF appears to be perfectly valid for heat conduction. 




10000 15000 20000 
t 

FIG. 3: Evolution of the local energies of a weakly coupled 
system as depicted in Fig. The initial state features one 
excited subsystem. Displayed are the predictions from the 
Kubo formula (solid lines) and the complete exact solution 
of the Schrodinger equation (points). The figure indicates 
diffusive transport in accord with the Kubo formula. 
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FIG. 4: Conductivity over frequency as calculated from the 
Kubo formula for a strongly coupled system as depicted in 
Fig.ED 



But are the criteria for diffusive transport from Sect.lVl 
fulfilled? An estimate for the correlation time is given 
by t c w 2ir/Se « 10 2 . From evaluating l|28|l one finds 
Td w 7 ■ 10 3 . Thus t c <C Td obviously holds. Furthermore 
2it/5uj w 6 • 10 3 , hence Td ~ 2ir/5u>, i.e., the criteria for 
diffusive transport are fulfilled. 

If, however, the interaction strength is such that the 
criteria from Sect. [V] are not fulfilled the transport be- 
havior ceases to be diffusive, i.e., the good agreement 
between the solution of (JSJ and the solution of the 
Schrodinger eq. vanishes. For example, for a model like 
the above one but with A = 5 • 10~ 4 one has Td ~ 70 and 
t c -C Td is not fulfilled. Fig. [S] shows the significant devi- 
ations of the evolution of the local energies from normal 
diffusive behavior as described by (JSJ . But the graph for 
k(u>) as calculated from the Kubo formula (see Fig. 
does not look essentially different from the above regular 
case. In particular there is no pronounced singularity at 
uj = as expected for the case of ballistic transport. Thus 
it is not obvious how the general transport behavior is to 
be found from simply evaluating the Kubo formula since 
there is no Td to be checked. The same is found for ex- 
tremely weak interactions. For example, for a model like 



FIG. 5: Evolution of the local energies of a strongly coupled 
system as depicted in Fig. Q The deviation of the points 
(Schrodinger equation) from the solid lines (Fourier's law, 
Kubo formula) indicate the breakdown of diffusive behavior 
and the validity of the Kubo Formula. 



FIG. 8: Gapless model to analyze diffusive transport: Two 
coupled subunits featuring n levels which are uniformly dis- 
tributed within the energy interval AE. 

0.04 | , , , , , , , 1 



se 

O 



-0.04 -0.02 0.02 0.04 



FIG. 6: Conductivity over frequency as calculated from the 
Kubo formula for an extremely weak coupled system as de- 
picted in Fig. Q 



the above one but with A = 10~ 5 one has w 1.8 • 10 5 
and Td «< 2it/5uj is not fulfilled. Consequently, the reg- 
ular transport behavior breaks down (see Fig. 0). But 
again, the graph for k(lo) as calculated from the Kubo 
formula does not look essentially different from the reg- 
ular case (see Fig. (BJ. 
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FIG. 9: Fourier transform of the current auto-correlation 
function for the model depicted in Fig. |H| (parameters see 
text). 



The other model is meant to show that the gapped 
spectrum of the subunits is dispensable for demonstrat- 
ing diffusive transport. It consists of subunits featuring n 
eigenstates distributed uniformly within an energy inter- 
val AE as depicted in Fig. |H| Here the interaction is cho- 
sen to be a complex random matrix on the full system's 
space without any restriction to a subspace. Neverthe- 
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FIG. 7: Evolution of the local energies of an extremely weakly 
coupled system as depicted in Fig. Q The deviation of the 
points (Schrodinger equation) from the solid lines (Fourier's 
law, Kubo formula) indicate the breakdown of diffusive be- 
havior and the validity of the Kubo Formula. 



FIG. 10: Evolution of the local energies for the model depicted 
in Fig. |H for an initial state featuring Ti > AE, T 2 < AE. 
Solid lines refer to the HAM prediction (Sj, points to the 
Schrodinger equation. The model clearly shows diffusive 
transport. 
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less, the interaction is supposed to be given by V = Xv 
with Tr{v 2 }/n 2 = 1. To keep the problem numerically 
manageable we consider only two subunits. We specify 
our model concretely by N — 2, n — 60, A = 5 • 10~ 3 , 
AE = 7. For this case we consider a ("pseudo-thermal") 
pure product initial state. The amplitudes of this state 
are chosen such that their magnitude squares obey Boltz- 
mann distributions with T\ = 40, T2 = 1. The phases 
of the amplitudes are chosen at random. This state has 
been chosen since for states with lower temperature dif- 
ferences it is hard to distinguish relaxation behavior from 
fluctuations. Since for this model various energy sub- 
spaces 77 yield different transport types and coefficients 
Dirf), the transport behavior for this initial state cannot 
simply be determined by directly applying the Kubo for- 
mula. It may, however, be analyzed within the framework 
described in Sect. I VII Since A 77 eventually determines the 
correlation time r c corresponding to the energy subspace 
77, one must, as eventually turns out, divide the energy 
scheme of this model in just two subspaces. Thus P(l) 
projects out all states with total energy E lower than AE, 
< E < AE, and P(2) the states with AE < E < 2AE. 
Based on those definitions one finds for the above initial 
state P(l) ~ 1, -P(2) « 0. Thus the relaxation behavior is 
controlled by D(l) as given by Fig.|5]shows the cor- 
responding C(l, uj) and we find (7(1, 0) = 0.036. Accord- 
ing to the above definitions we have t c « 2tt/AE ks 0.9. 
From numerically evaluating we find « 40 and thus 
the corresponding condition for normal energy diffusion 
t c < r d is fulfilled. With 5lj = 0.3 we get 2-k/Slj « 21 
which is of the same order of magnitude as , which also 
indicates diffusive transport. From numerical evaluation 
we furthermore find Tr{po(l)S 2 } = 2.03. Those numbers 
eventually yield D(l) — 0.0139. Indeed, as displayed in 
Fig. ^3 there is good agreement of the solution of (J3J) 
based on this D(l) with the full numerical solution of 
the corresponding Schrodinger equation. 



as weakly coupled identical many- level subunits. We ar- 
gued that this description may apply to a large class of 
systems. For simplicity, we concretely analyzed modular 
chains and rings with next neighbor couplings. With- 
out making any reference to external forces or exploiting 
the hypothesis of local equilibrium we showed that those 
systems may or may not exhibit normal heat transport, 
but if they do, the conductivity is correctly described by 
the Kubo formula. This is in accord with Refs>iiiii^ 
where the Kubo Formula has been evaluated for con- 
crete systems and the results have been counter-checked 
by either experiments or other theoretical methods. This 
is furthermore in accord with R e f Sl 22 4 22iiffii2ii22i where it 
is shown that one-dimensional systems may or may not 
exhibit diffusive heat transport. 

We also suggested general criteria to decide whether 
or not such modular systems exhibit normal transport. 
Those criteria are established on the basis of the concrete 
form of the subunits and their mutual interactions. We 
found, however, that the question cannot be decided only 
by evaluating the Kubo Formula. To check our theoreti- 
cal results we introduced some examples for concrete, fi- 
nite, modular, chainlike systems. For those examples we 
numerically solved the time-dependent Schrodinger equa- 
tion which yields the energy transport dynamics. Those 
dynamics are in accord with our above mentioned theo- 
retical results. 



The results support the view that thermodynamic be- 
havior might, under specific conditions on the system, 
emerge directly from quantum mechanics^22i24. These 
conditions do not necessarily include a many particle 
limit. 



VIII. SUMMARY AND CONCLUSION 

We investigated the energy transport behavior of mod- 
ular quantum systems, i.e., systems that can be described 
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